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Polymers with both soluble and insoluble blocks typically self-assemble into micelles, aggregates 
of a finite number of polymers where the soluble blocks shield the insoluble ones from contact with 
the solvent. Upon increasing concentration, these micelles often form gels that exhibit crystalline 
order in many systems. In this paper, we present a study of both the dynamics and the equilibrium 
properties of micellar crystals of triblock polymers using molecular dynamics simulations. Our 
results show that equilibration of single micelle degrees of freedom and crystal formation occurs 
by polymer transfer between micelles, a process that is described by transition state theory. Near 
the disorder (or melting) transition, bcc lattices are favored for all triblocks studied. Lattices with 
fee ordering are also found, but only at lower kinetic temperatures and for triblocks with short 
hydrophilic blocks. Our results lead to a number of theoretical considerations and suggest a range 
of implications to experimental systems with a particular emphasis on Pluronic polymers. 



INTRODUCTION 

Micelles of multi-block polymers are finite aggregates, 
typically around fifty polymers or less, where the insol- 
uble blocks shield the soluble ones from contact with 
the surrounding solvent. Depending on control variables 
(temperature, polymer concentration, pH, etc..) micelles 
may self-assemble into gels that exhibit long-range or- 
der such as bcc, fee, hep or other, more unusual, crys- 
tals. Micellar crystals exhibit a number of unique prop- 
erties that have made them extremely attractive for fun- 
damental studies as well as for applications [1, 2]. Ex- 
tensively studied experimental systems include aqueous 
solutions of Pluronics (also known as Poloxamers), ABA 
triblocks where A is Polyethylene oxide (PEO) and B is 
Polypropylene oxide (PPO) [3, 4, 5] and inverted Pluron- 
ics, where the A blocks are PPO and the central block 
B is PEO [5, 6] as well as non-aqueous systems such as 
Polystyrene-Polyisoprene (PS-PI) diblocks in decane [7] 
and other solvents [9, 10, 11, 12]. 

Micelles in solution are highly dynamical entities 
with polymers continually being absorbed and released 
through time. Therefore, a micellar crystal has a con- 
siderably intricate structure, where the long range order 
remains stable as the individual polymers are constantly 
hopping from one micelle to the next. Theoretical ap- 
proaches such as density functional or mean field theory 
[2, 8, 13, 14, 15, 16, 17] directly study the ordered mi- 
celles and ignore the dynamical degrees of freedom of the 
polymers. Studies using molecular dynamics (MD) have 
the advantage of providing a reasonably realistic descrip- 
tion of the dynamics, thus allowing the investigation of 
the role of single polymer degrees of freedom. In contrast 
with other approaches, MD also offers the important ad- 
vantage that no assumptions need to be made about what 
is the possible thermodynamic state of the system. 



The goal of this study is to predict phase diagrams of 
triblock polymers using MD simulations and to gain an 
understanding of the dynamics of micellar crystal forma- 
tion. Because of our ongoing interest in Pluronic sys- 
tems in aqueous solutions [18], we examine systems of 
A n B m A n triblocks, where the A beads are hydrophilic 
and B beads are hydrophobic. Although there have been 
a number of previous studies of multiblock copolymers 
in solution using MD, see Ref. [19] for a recent review, 
the prediction of crystalline structures presents substan- 
tial difficulties. Experimentally, it is well known that the 
approach to thermodynamic equilibrium is slow in these 
systems, with time scales of the order of minutes or hours. 
Therefore, even with suitably coarse-grained models, the 
long time scales involved provide a considerable challenge 
for MD simulation studies. 

In this paper, we provide a detailed investigation of 
self assembled micellar crystals using MD. We aim to 
understand the mechanism by which they form, predict 
their range of stability and elucidate their static and dy- 
namic properties. We also present an in-depth study on 
the challenges associated in reaching the thermodynam- 
ically stable state using MD and a successful strategy to 
overcome them. 

MODEL AND SIMULATION DETAILS 
Simulation Details 

Systems of polymers are modeled by coarse-grained 
beads in an implicit solvent. Ref. [18] provides a more de- 
tailed justification than the outline provided here. Indi- 
vidual polymers are A n B m A n symmetric triblocks, where 
A beads are hydrophilic and B beads are hydrophobic. 
All systems in this work are monodisperse with fixed val- 
ues for n and m. Non-bonded pair potentials consist 
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of a standard attractive Lennard- Jones potential for hy- 
drophobic interactions 
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and a purely repulsive potential for hydrophilic interac- 
tions 
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Pair potentials are cutoff to zero at r = 3.0a. All beads 
have the same mass M. M, e and a are arbitrary and 
uniquely define the units of all numbers in the simula- 
tions. Neighboring atoms in the polymer chain are con- 
nected together with a simple harmonic potential. The 
packing fraction <f>p of the polymers is given by 
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where N po \ y is the number of polymers in a box of lin- 
ear dimensions L and N mon — 2n + m is the number of 
beads in each polymer. Simulation boxes are cubic with 
periodic boundary conditions. For this work, we use a 
time step of At — 0.005-^/ ma 2 je. All simulations were 
performed using the LAMMPS software package [20] in 
the NVT ensemble via Nose-Hoover dynamics [21]. 

The temperature sensitivity of Pluronic systems is a 
reflection of the underlying strong temperature depen- 
dence of the hydrophobic effect. To model temperature 
dependent phases in such systems, the solvent quality for 
the A beads must change [18]. In this paper, we keep the 
solvent quality fixed and vary only the kinetic tempera- 
ture when needed to equilibrate the systems properly. 

Recently, a different implicit solvent model has been 
developed for describing Pluronic systems [22], where 
the pair potentials for the coarse-grained simulation are 
fitted to results obtained from all-atom MD and quan- 
tum chemistry simulations. Although the coarse-graining 
is different (each monomer A,B represents one PEO or 
PPO monomer) than in our model, the resulting poten- 
tials are quite similar to the ones in this work, Equations 
1 and 2, with the only significant difference that the val- 
ues of a are different between the two A and B monomers, 
and Uaa shows a minor maximum. 



Observables 

A number of observable quantities are monitored for 
every recorded time step during a simulation run. The 
micelles themselves arc identified by the same algorithm 
used previously in Rcf. [18]. Any hydrophobic beads 
within a distance of r cut = 1.2c of one another are identi- 
fied as belonging to the same micelle. Identified micelles 
containing less than 3 polymers are typically free poly- 
mers in the process of being transferred from one micelle 



to another and are removed from further consideration. 
Observables such as statistics of micelle aggregation num- 
ber, gyration tensor, center of mass, and micelle lifetime 
are calculated and examined for every simulation per- 
formed in this work. Methods used to calculate these are 
described in Ref. [18]. 

The structure factor S(q) is calculated over the center 
of mass coordinates fi of all N m i c micelles in the system 
using the formula 



S(q) = C ( 
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with the components of q as multiples of ^ due to the 
use of periodic boundary conditions. The peaks in S(q) 
are then used to reconstruct the full 3D real space lat- 
tice basis, if it exists. In this manner, S(q) is not be- 
ing used to simulate real scattering intensities as may be 
obtained by X-ray experiments, but as a mathematical 
order parameter to discriminate between the different or- 
dered structures that may be present. For convenience, 
the normalization Cq is chosen so that S(q = 0) = 1. A 
more sophisticated treatment that allows continuum val- 
ues of q, suitable for quantitative comparisons with X-ray 
experiments, has been recently introduced [23], but we do 
not use it here. 

Individual polymers are constantly hopping from one 
micelle to another. This is quantified over the entire sim- 
ulation box as an overall rate of polymer transfer, rpx, by 
examining contiguous simulation snapshots. Sets of in- 
dexed polymers belonging to each micelle are compared 
between the snapshots to find the number of polymers 
transferred. The rate rpx is then expressed as a fraction 
of the number of polymers in the box transferred per one 
million time steps. A polymer that is transferred out and 
back to the same micelle between snapshots will not be 
counted by this analysis, so snapshots are recorded ev- 
ery 100,000 time steps to minimize undcrcounting. At 
UbT/e — 1.2, a typical micelle only loses/gains one poly- 
mer per ten snapshots recorded. 

Radial distributions of the beads surrounding micelles 
are also of interest. These are calculated by creating a 
histogram with bin width dr and then counting the num- 
ber of beads belonging to a micelle N count (r) that fall 
between r and r + dr, where r is the distance of the bead 
from the center of mass of the micelle. Good statistics 
require averaging this histogram over all micelles in the 
simulation and over all time steps after the micellar crys- 
tal has formed. The average histogram is transformed 
into a radial density distribution of beads around a mi- 
celle by calculating the packing fraction of beads in each 
bin 
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where i refers to either A or B beads. We use dr — 
0.2a to balance smooth graphs with the need for long 
simulation runs to obtain detailed statistics. 

MICELLAR CRYSTALS STUDIED USING MD 

There are two major challenges faced in using MD to 
determine equilibrium phases. First, the simulation must 
last longer than any of the relaxation times in the system. 
Second, the simulation box size must be chosen properly 
to avoid finite size effects. The first problem is related 
to the kinetic temperature at which the system is run. 
The second problem can become particularly severe for 
simulating crystals with three dimensional order, where 
the incorrect choice of an even large box size L can force 
the system into very distorted ordered phases. These two 
issues are addressed systematically using simulations of 
the A10-B7A10 polymer. It provides a coarse-grained de- 
scription of one of the most extensively studied Pluronics, 
F127 [3]. The conclusions of this study lead to a general 
methodology valid for any other polymer. 

140 1 . . 1 . 1 




time step /10 1 



FIG. 1: Mean squared displacement averaged over all beads 
in the simulation box. This data was obtained from the initial 
simulation runs of the A10B7A10 polymer at iVpoiy = 500 
and 4>p = 0.20. The origin of the time axis is 30 million 
time steps which indicates that the recording of these results 
began after the system reached equilibrium. In the case of 
ksT/e — 1.0 this equilibrium is a metastable state, and the 
beads are not diffusing. The simulation run performed at 
knT/e — 1.2 formed a fee lattice around time step 10 million 
which persisted until the end of the run at 35 million, and the 
mean squared displacement shows a characteristic diffusive 
behavior. 



Micelle crystallization and kinetic temperature 

Initial simulations are performed with a system size of 
IVpoiy = 500 at a kinetic temperature = 1.0 and 




FIG. 2: Snapshot of a A10B7A10 polymer simulation run at 
4> P = 0.20, k B T/e = 1.2, 7V po i y = 500, taken after the fee 
lattice formed. A beads are represented by orange spheres, 
and are shown with a reduced radius so they do not obscure 
important details. Orange lines indicate bonds between these 
beads. B beads are shown in blue with a radius of 0.6a. 
Large yellow spheres are placed on the lattice reconstructed 
from S(q). Every yellow sphere is sitting on a micelle, visually 
confirming a perfect fee crystal. The A beads have been re- 
moved around a single unit cell of the lattice and yellow lines 
added to guide the eye. All snapshots are generated using 
PyMol [24]. 



with concentrations <j) P of 0.05, 0.10, 0.15, 0.20 and 0.25 
run over 30 million time steps each. In all cases, the 
initially randomly placed A10B7A1Q polymers aggregate 
into micelles quickly in a few thousand time steps. 

Visual examination indicates that concentrations at 
and above <fip = 0.15 are strong candidates for micel- 
lar crystals. Micelles are locked in place and only move 
a few a from their average positions. The analysis of the 
order parameter S(q), however, indicates no long-range 
ordered structures exist in any of these initial simula- 
tions. An inspection of the polymer transfer shows that 
it remains negligible throughout all simulations. This 
is confirmed by calculating the mean squared displace- 
ment ((r(t) — r(0)) 2 ) averaged over all beads in the sys- 
tem. Figure 1 clearly shows non-diffusive behavior at 
ksT/e = 1.0. These results suggest that the individual 
micelles are quickly frozen in a configuration that is not 
representative of equilibrium, thus preventing the entire 
system from reaching thermal equilibrium. 

The lack of equilibration of micellar degrees of freedom 
suggests subsequent runs at a larger kinetic temperature 
ksT/e — 1.2. A single simulation run at (ftp = 0.20 
formed a textbook fee lattice after about 10 million steps, 
and is shown in Figure 2. Throughout the duration of the 
run, the rate of polymer transfer was substantial, about 
7% of the polymers in the box every million time steps, 
even after equilibrium is reached. Playing a movie of the 
simulation shows that after the lattice formed, micelles 



4 




126 



10 20 30 

time step /10 6 



40 



FIG. 3: Examples of the behavior of iV m i c during a simulation 
run of the AioBrAio polymer. Results are shown here for two 
independent runs with <j>p = 0.20, kpT/e = 1.2 and N po i y = 
2134. The origin of the time axis (t=0) indicates the time 
step where the entire simulation started from a random initial 
configuration. Note how both systems eventually plateau at 
the same state N m ic = 128, though it occurs at different times. 



do not appear to move except by vibrating about their 
average positions. However, while the micelles appear 
static, polymers are constantly being exchanged among 
the micelles, so that any individual polymer will even- 
tually explore the entire simulation box. This is inde- 
pendently confirmed by the analysis of the mean squared 
displacement of beads in Figure 1, which shows a classic 
diffusion result at ksT/e — 1.2. 

The behavior of the number of micelles N m - lc in the box 
as a function of simulation time is of particular interest. 
In simulation runs where micellar crystals are found, the 
system reaches a plateau where N m - lc remains constant, 
see Figure 3 for an example. Despite their dynamic char- 
acter even at equilibrium, -ZV m i C then remains constant 
for the duration of the run. This correlation is a general 
feature in all simulation runs performed. Every single 
one that leads to a stable plateau in iV m i c as a function 
of time formed a micellar crystal confirmed by peaks in 
the order parameter S(q). 



Box size effects 

The influence of the simulation box size choice is as- 
sessed by running additional simulations with N po i y = 
500, 600, 800, 900, and 1000 with fixed concentration 
cf)p = 0.20 and — 1.2. Several different random ini- 
tial configurations are used at each system size to ensure 
repeatability. Table I summarizes the results of all these 
simulation runs. 

Interestingly, the additional three simulation runs at 
•Wpoly = 500 do not exhibit fee structures. Instead, each 




FIG. 4: Snapshot of a A\§B-jA\§ polymer simulation run at 
4>p = 0.20, k B T/e = 1.2, and iV po i y = 600 taken after the 
lattice formed. The resulting lattice in this system is a body- 
centered tetragonal with c/a — 1.5. Coloring conventions are 
identical to Figure 2 



No. configs 


-/Vpoly 


iVmic 


Ordering 


<^a gg ) 


1 


500 


32 


textbook fee 


15.625 


3 


500 


unstable 


none 




1 


600 


36 


none 


16.67 


3 


600 


36 


distorted bec 


16.67 


5 


800 


48 


distorted bec 


16.67 


5 


900 


54 


textbook bec 


16.67 


1 


900 


55 


distorted bec 


16.36 


2 


1000 


unstable 


none 




4 


1000 


60 


distorted bec 


16.67 



TABLE I: Summary of results obtained from initial test runs. 



of them remains unstable with iV m ; c never achieving a 
constant plateau and S(q) is devoid of peaks. Larger 
system sizes do reach stable ordered structures with con- 
stant N m i c and numerous peaks in S(q). Most of the 
structures that occur appear bec when examined visually, 
but a detailed analysis of the lattices indicates that many 
of them are distorted. One obviously distorted structure 
is depicted in Figure 4 where the lattice is body centered 
tetragonal with c/a = 1.5. 

Other distortions include body centered tetragonal 
with c/a — 1.054 for N po \ y = 800 and a lattice that 
appears to be almost exactly bec when N po \ y = 1000, 
except that the central micelle in the unit cell is shifted 
slightly from the true center. Lastly, a textbook bec lat- 
tice is formed in the simulation runs with N po i y = 900. 

Examining the average aggregation number leads to a 
very illuminating result. It is found that for almost all 
simulations these numbers are identical. This indicates 
that the average micelle aggregation number is indepen- 
dent of the box length (at a fixed 4>p) & n d is only a func- 
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No. configs 


Targeted 


A^poly 


N ■ 

1 v mic 


Ordering 


1 


16 micelle bcc 


267 


16 


textbook bcc 


4 


108 micelle fee 


1800 


unstable 


none 


4 


128 micelle bcc 


2134 


128 


textbook bcc 


3 


250 micelle bcc 


4168 


unstable 


none 



TABLE II: Summary of simulation results from testing the 
algorithm on A10B7A10. 



tion of the polymer structure, kinetic temperature and 
concentration. 

Assuming without proof that either the fee or the bcc 
lattice represents the real thermodynamic equilibrium of 
the system, the previous observation then suggests a way 
to generate magic numbers of polymers to simulate equi- 
librium states free of finite size effects. In a bcc lattice, 
each cubic unit cell contains C — 2 micelles, while the 
fee lattice contains C = 4. Therefore, in order to obtain 
a bcc or fee lattice with M by M by M unit cells in a 
cubic simulation box, the number of polymers needed to 
achieve this is given by 



iV poly = CM 3 (iV agg ). 



(5) 



MD Simulations without finite size effects 



An algorithm to simulate micellar crystals without fi- 
nite size effects follows very naturally from the previous 
results, and is summarized in the following steps. 

1. Concentration selection: The concentration must 
be chosen high enough so that micelles pack into a 
potential crystalline state. 

2. Temperature selection: The temperature should 
be chosen large enough to ensure a significant rate 
of polymer transfer, but low enough that the micel- 
lar crystals are not in a disordered phase. As a rule 
of thumb, we have been using a polymer transfer 
around 10% of the polymers in the box every mil- 
lion steps. 

3. System size selection: Calculate the average mi- 
celle number via test simulations and use Equa- 
tion 5 to determine the final system sizes to perform 
simulations on. 

4. Ensure reproducibility: The formation of micellar 
crystals is a stochastic process so several simula- 
tions with different initial configurations must be 
run. 

The advantage of this algorithm is that steps 1-4 can be 
accomplished with relatively modest computer resources 
on small system sizes, leaving the production runs with 
large polymer numbers as the only computationally in- 
tensive calculations. 




FIG. 5: Snapshot of a A10B7A10 simulation run at <f)p = 0.20, 
ksT/e = 1.2, and iVpoiy = 2134, taken after the bcc lattice 
formed. Coloring conventions are identical to Figure 2 



Micellar crystals of general A n B m A n triblocks 
AiqBtAio 

Further simulation runs are carried out on the 
A W B 7 A 10 polymer system to test the algorithm. These 
additional simulations target 16, 128 and 250 micelle bcc 
configurations, along with 32 and 108 micelle fee ones. 
The simulations are summarized in Table II. None of the 
simulation runs targeting fee phases ever reach a stable 
number of micelles and correspondingly, S(q) indicates 
there is no long-range order. On the other hand, both 
the 16 and 128 micelle bcc configurations are perfect and 
completely reproducible with the expected ordering con- 
firmed by the structure factor. Figure 5 shows a snapshot 
of the 128 micelle bcc phase after it forms. Further dis- 
cussion on the calculated structure factors are presented 
below. Larger simulations attempting the formation of 
a 250 micelle bcc crystal never resulted in a stable iVmi 
plateau. 

-4.20-Bi4j4.20 

The A20B14A20 polymer is also studied as it provides 
a closer representation to the real Pluronic F127, since 
there are twice as many monomers in a polymer and the 
level of coarse-graining is less, as discussed in Ref. [14]. 
Applying the algorithm developed above to search for 
micellar crystals, the concentration was selected at <fip — 
0.20 and the kinetic temperature at ksT/e = 1.9. Ini- 
tial simulation runs establish that (N a g g ) — 22.5. The 
results, summarized in Table III, are similar to those for 
the smaller AiqB^Aiq polymer, with the bcc structure 
being the most commonly found. The fee phase again 
only appears in a single simulation run and is not repro- 
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No. configs 


Targeted 


Apoly 


N mic 


Ordering 


No. configs 


Targeted 


Apoly 




Ordering 


2 


16 micelle bcc 


360 


16 


textbook bcc 


4 


32 micelle fee 


444 


32 


textbook fee 


1 


32 micelle fee 


720 


32 


textbook fee 


1 


32 micelle fee 


444 


unstable 


none 


1 


32 micelle fee 


720 


35 


distorted bcc 


2 


108 micelle fee 


1500 


unstable 


none 


2 


32 micelle fee 


720 


unstable 


none 


2 


108 micelle fee 


1500 


103 


distorted bcc 


1 


54 micelle bcc 


1215 


54 


textbook bcc 


1 


108 micelle fee 


1500 


104 


distorted bcc 


3 


54 micelle bcc 


1215 


55 


distorted bcc 


TABLE V: Summary of simulation results from testing th 
algorithm on A^BjA^ after cooling to fcsT/e = 1.07. 


2 


108 micelle fee 


2430 


unstable 


none 



TABLE III: Summary of simulation results from testing the 
algorithm on A2qB\4,A2o at ksT/e = 1.9. 



DYNAMICS OF CRYSTAL FORMATION 



No. configs 


Targeted 


A r poly 


Amic 


Ordering 


5 


16 micelle bcc 


222 


16 


textbook bcc 


5 


32 micelle fee 


444 


unstable 


none 


5 


54 micelle bcc 


750 


54 


textbook bcc 


5 


108 micelle fee 


1500 


unstable 


none 


3 


128 micelle bcc 


1778 


128 


textbook bcc 


1 


128 micelle bcc 


1778 


136 


distorted bcc 


3 


250 micelle bcc 


3472 


unstable 


none 


4 


432 micelle bcc 


6035 


unstable 


none 



TABLE IV: Summary of simulation results from testing the 
algorithm on A§B~iA(, at ksT/e = 1.1. 



Molecular dynamics simulations not only allow the pre- 
diction of equilibrium phase diagram, but also describe 
the dynamics of micelles as they evolve towards thermal 
equilibrium. In the simulations presented previously, mi- 
celles quickly form from a completely random configura- 
tion of polymers, and then after a long simulation time 
(10 to 20 million time steps) order into a micellar crystal. 
We now present a quantitative picture of the dynamics 
of the polymers and micelles as the system approaches 
equilibrium. All results in this section are obtained from 
the analysis of the simulations of the Ai^B^A^ polymer. 

Polymer transfer is an activated process 



duciblc. 



A 6 B 7 A 6 

Given the prevalence of bcc lattices so far, a poly- 
mer with shorter ^4-blocks {AqB^Aq), expected to form 
more crew-cut micelles and hence more prone to as- 
semble into a fee lattice, is also investigated. Apply- 
ing the algorithm developed above, the concentration 
is selected at </>p = 0.15, the kinetic temperature at 
ksT/s = 1.1, and the average aggregation number is 
found to be (N agg ) — 13.875. Again, the bcc micellar 
crystal is most commonly found, as shown in the results 
in Table IV. No fee phases are found at this kinetic tem- 
perature. 

However, the rate of polymer transfer at fc^T/e =1.1 
is larger than 10% per million steps, so some addi- 
tional simulation runs are performed at a slightly reduced 
temperature ksT/e — 1.07, where polymer transfer is 
slightly reduced. In this case a completely reproducible 
fee structure is found for relatively small system sizes, 
but no fee lattices could be stabilized for larger ones as 
shown in the results in Table V. 



Polymer transfer plays an important role in achiev- 
ing the formation of micellar crystals in the simulations 
discussed above. If there is too little, the single micelle 
degrees of freedom do not reach equilibrium. Too much 
pushes the system into a disordered state. Moreover, the 
rate of polymer transfer is extremely sensitive to the ki- 
netic temperature. An equilibrated A^iy = 2134 bcc 
micellar crystal was taken as an initial configuration for 
additional simulations that continued with fc^T/e rang- 
ing from 1.0 to 1.3. Figure 6 shows the results. The rate 
of polymer transfer rpx starts near at ksT/e = 1.0 and 
increases exponentially, following an Arrhenius form 

/ AG K 

r PT = r exp(- — — ) . (6) 
k B T 

That is to say that polymer transfer is an activated pro- 
cess. Curve fitting, we find 

AG* « 10e . (7) 

Figure 6 also includes results from simulations per- 
formed using a Langevin thermostat [25] . It controls the 
temperature by adding an additional force to every par- 
ticle F — — ■yv + F ran d where the magnitude of the ran- 
dom force -Frand and 7 set the temperature through the 
fluctuation dissipation theorem [25]. The results for two 
different values of 7, which are plotted in Figure 6, show 
the dependence of the polymer transfer for different drag 
coefficients. 
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k B T/e 

FIG. 6: Polymer transfer rpx versus temperature calculated 
from simulation runs of the A10B7A10 polymer at <f)p = 0.20 
and iVpoiy = 2134. All simulations start from an already equi- 
librated bcc phase. Results are included for the Nose-Hoover 
thermostat and Langevin thermostat with two different val- 
ues of 7. The inset plots the same data as a plot of log 10 (rpT) 
versus e/figT) to show that the slopes of the resulting lines 
(—AG*) are universal. 



In Figure 6, the slope of the lines in the inset plot are 
universal for both thermostats and both values of 7. The 
universality of this value implies that the calculated AG' 
is the free energy barrier between a polymer in a micelle 
to a transition state between micelles. While their slopes 
are universal, the y-intercepts (related to ro) should de- 
pend on the diffusion coefficient of the hydrophilic beads, 
which, from the Einstein relation [26] is inversely pro- 
portional to the drag coefficient. This is reflected in the 
offsets of the various plots in Figure 6, which are clearly 
different . 

The difference in free energy between a polymer within 
a micelle and in the transition state entirely surrounded 
by solvent and hydrophilic beads can be roughly esti- 



mated as AG" 



L exp 



e, where 



'•cxp 



is the number 



of hydrophobic beads exposed to solvent. If the length 
of hydrophobic block is low, m 



exp 



n. Thus for the 
AiqBiAiq polymer analyzed here, it is expected that 
AG" ~ 7s, which is consistent with the measured value. 
Hydrophobic beads in polymers with a longer hydropho- 
bic block are expected to form a globule in the transi- 
tion state, leading to a slower increase of AG" going as 
Assuming that ro remains constant as 



m. 



in 



2/3 



m increases, then maintaining the same rate of polymer 
transfer rpx will require increasing the kinetic temper- 
ature by the same factor. These considerations agree 
remarkably with the kinetic temperature of fc^T/e = 1.9 
selected for the A20B14A20 polymer simulations, 

( rnaeVL ) ^ • k B T wcv /e = ( ~ • 1.2 = 1.9 

Wprcv / V 7 / 
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FIG. 7: Number of micelles in the ordered phase N OI d as a 
function of time for a single simulation run of the A10B7A10 
polymer at <f)p — 0.20 and temperature ksT/e = 1.2. The 
origin of the time axis (t=0) on this plot indicates the time 
step where the simulation was started from a random config- 
uration. 



However, as the size of the hydrophobic block grows, 
the transition state should also have an additional con- 
tribution due to the free energy cost of passing a globule 
through the brush formed by the hydrophilic coronas. So 
this simple estimate should eventually break down. 

Implicit in our arguments is that polymer transfer ac- 
counts for the diffusive behavior in Figure 1. Therefore, 
the diffusion coefficient can be roughly estimated from 
the time tpt it takes a polymer to travel the nearest 
neighbor micelle distance a^. Comparing this estimate 
with the fit to Figure 1, a good agreement is obtained: 



((r(t)-r(O)) 2 



1.3 • 10" 



Tpt 



a L r PT 



= 1.1 • 10 



-5 1 



At 
At 



(8) 



Dynamics of micellar crystal formation 

The structure factor S(q) is sufficient as an order pa- 
rameter to determine if the entire system is in an ordered 
state, but it reveals no information about the dynamics 
before that final state is formed. For this, we turn to the 
bond order analysis [27] and apply it to all the micelles 
at every time step in simulation runs performed on the 
AiqBtAio polymer with iVpoiy = 2134. In short, at any 
given time step where a micellar crystal may have par- 
tially formed, the bond order analysis identifies those mi- 
celles belonging to the ordered and the disordered phases. 

One simple way to examine the results is to count the 
number of ordered micelles AT or( j in the simulation box 
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at every time step. An example from one simulation 
run is shown in Figure 7, plots for all other simulation 
runs performed are qualitatively very similar, though the 
ordered phase appears at different times. 

After starting the simulation shown in Figure 7 from 
a random configuration, micelles quickly form in only a 
few thousand time steps and single micelle degrees of free- 
dom are equilibrated before one million time steps have 
passed. The system then explores configuration space 
as polymers are transferred and micelles travel anywhere 
from a few a up to lOOcx without any micelles appear- 
ing ordered until time step 10 million. At this point, 
the number of ordered micelles grows very quickly over 
the next one million steps until all 128 micelles are in 
the bcc lattice and remain there for the duration of the 
simulation. 

During this short time span of iV or d growth, only 7% of 
the polymers in the box are transferred between micelles, 
such a small amount that it cannot fully account for the 
ordering of all the micelles. During the same time in- 
terval, a detailed analysis shows that some micelles move 
significantly (up to 10a) before the ordered phase finishes 
forming. This may suggest that micellar crystal forma- 
tion is a two step process, where first individual micelles 
are equilibrated by polymer transfer followed by a second 
step where polymer transfer becomes irrelevant and the 
actual crystal grows via the movement of micelles. 

However, some additional simulations performed to 
test this hypothesis do not support it. Single micelle 
degrees of freedom in these tests are first equilibrated at 
UbT /e = 1.2 for 5 million time steps and then the ki- 
netic temperature is quenched to k^T/e = 1.0 to signif- 
icantly reduce polymer transfer and allow micelle move- 
ment. None of these simulation runs resulted in the for- 
mation of an ordered phase. We therefore conclude that 
a significant amount of polymer transfer remains a crit- 
ical component in the actual growth of ordered micellar 
crystals. 
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FIG. 8: Structure factor calculated after the lattice formed 
for the A10-B7A10 polymer simulation run at 4>p = 0.20, 
k B T/e = 1.2, and iV po i y = 2134. The full 3D S(q) is plot- 
ted as a scatter plot of S(q) versus \q\. The multiplicity of 
the various peaks can be seen. Vertical dotted lines indicate 
the location of identified peaks, and their positions relative 
to q* = 4 • 27r/L are also noted (the factor of 4 is included 
because there are 4 unit cells along the box length L). 
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PROPERTIES OF MICELLAR CRYSTALS 
Structure factor 

We calculate the structure factor (Equation 4) for the 
equilibrium state of every simulation run performed and 
use it as an order parameter to determine if a micellar 
crystal is present. In those systems where the crystal does 
form it is a perfect single crystal and, correspondingly, a 
large number of peaks can be identified in S(q) as shown 
in Figure 8. Peaks occur in reciprocal space at discrete 
points q — G, where G are the reciprocal vectors of the 
corresponding lattice. 

Peaks evident in Figure 8 decrease in magnitude for 
larger values of G. This damping is expected to follow a 



FIG. 9: Fluctuations in micelle positions vs. temperature 
calculated for a A\qBtAu) polymer simulation run at <j>p = 
0.20, k B T/e = 1.2, and 7V poly = 900. The micellar crystal 
was formed at k B T/e = 1.2 and then cooled to the target 
temperature without disrupting the lattice. Simulation runs 
heated to a higher temperature do disrupt the lattice and Ar 2 
becomes ill-defined. Fluctuations on the y-axis are plotted as 
a ratio relative to a B , the nearest neighbor distance in the 
lattice. 

Debye- Waller factor [28] approximately described as 

S(q = G)cxexp(-(Ar 2 )|G| 2 /3) , (9) 

where (Ar 2 ) is the mean square displacement of the mi- 
celle center of mass from its ideal lattice position. A 
curve fit, shown in Figure 8, is in excellent agreement 
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with Equation 9. The Lindemann ratio Jl is defined as 

v / (Ar2) = / i a i (10) 

where az, is the nearest neighbor distance between mi- 
celles. Using the parameters of the curve fit yields 
h « 0.14. 

The Lindemann parameter can alternatively be 
computed directly by measuring the mean square dis- 
placement of micelles about their average lattice posi- 
tions. The results, shown in Figure 9 agree remarkably 
well with the estimate from the Debye- Waller factor. Fur- 
thermore, it has been established empirically that ap- 
proximately at /i = 0.13 solids melt into disordered 
states [29] , a result that is also supported from our simu- 
lations as we observed no micellar crystals form at kinetic 
temperatures greater than ksT/e = 1.2 
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FIG. 10: Radial density distribution for two nearest neigh- 
bor micelles superimposed with a separation of the nearest 
neighbor distance in the bcc lattice. The y-axis plots the vol- 
ume fraction of the different beads belonging to a micelle in 
the local environment around the micelle. The results were 
calculated from a A10B7A10 simulation run at <f)p — 0.20, 
ksT/e — 1.2, iVpoiy = 2134 and averaged over all micelles 
and time steps after the micellar crystal has formed. 



Structure of the lattice of micelles 

Micelles in the AiqBtAiq polymer system are arranged 
with the centers of mass sitting on a bcc lattice with a 
nearest neighbor spacing of = 11.53c This num- 
ber remarkably agrees with the length of the polymer if 
stretched completely into a straight line across the diam- 
eter of a circle from the center of one nearest neighbor mi- 
celle to the opposite one, N nlon -r = 27-0. 83cr = 22 Ala ~ 
2dL, where r$ — 0.83cr is the equilibrium bond length. 
This suggests that polymers are maximally stretched, ei- 
ther along the diameter or in a bent configuration. De- 
tailed visual examinations of simulation snapshots indi- 



cate that there are a significant number of polymers in 
the bent configuration, but there are also some in the lin- 
ear extended configuration. The micellar cores are liquid- 
like, and over time, a given polymer constantly switches 
from linear to bent configurations. 

The previous arguments also suggest a significant 
amount of overlap of the coronas between two neighbor- 
ing micelles in the lattice. To examine this more quanti- 
tatively we calculate the micelle density distribution as a 
function of the radius averaged over all micelles and all 
time steps in the simulation after the micellar crystal has 
formed. Figure 10 shows the results, confirming the over- 
lap. Hydrophilic A beads from one micelle explore the 
solvent until occasionally bumping into the hydrophobic 
core of one of the nearest neighbor micelles. 

A remarkable aspect of micellar crystals found in this 
work is the stability of the average aggregation num- 
ber. If the polymers within the micelle are maximally 
stretched, the core of the micellar radius is R c = totoct/2, 
so the aggregation number can be estimated from 



(N agg ) = 



47r(mr /2) 3 (7 3 ir 



,2„2 



mro<7 



6 



(11) 



For the AioBjAiq this yields (N agg } = 17.7, in good 
agreement with the simulation results in Table I. Aggre- 
gation numbers for the other polymers simulated do not 
agree, implying that the polymers in those systems are 
not maximally stretched. 

CONCLUSIONS 
Summary of results 

Cubic micellar crystals of A n B m A n polymers form in 
MD simulations at sufficiently high concentrations. In 
order to form the crystals, high enough kinetic tempera- 
tures are needed to enable polymer transfer between mi- 
celles, which is critical for equilibrating the system. The 
polymer transfer process is activated and described by 
transition state theory. It results in an apparent diffu- 
sive behavior of the individual polymers while the lattice 
of micelles remains stable. Excessive polymer transfer 
at even higher kinetic temperatures triggers a disorder 
phase transition to a micelle liquid. 

In the process of forming the ordered phase, the sys- 
tem spends a long time in a micellar liquid phase equi- 
libration period where no ordered nucleates are present. 
Eventually, a large nucleation event takes place and the 
micellar crystal then grows very quickly, filling the entire 
simulation box in a relatively short time span. During 
this growth period, polymer transfer and movement of 
micelles are both crucial in the formation of the final 
micellar crystal. The preferred ordering near the disor- 
dered transition for all triblocks studied in this system is 
the bcc lattice. Only at lower kinetic temperatures and 
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FIG. 11: Summary of the phase diagram encompassing all 
simulated triblocks A n B m A n (all forming cubic phases with 
long range order). Near the disorder transition, bcc is al- 
ways favored and fee lattices only begin to appear at lower 
temperatures. At even lower temperatures polymer transfer 
becomes negligible and MD would require prohibitively long 
simulations to reach equilibrium. 

for polymers with short hydrophilic groups (AqB^Aq at 
UbT/e = 1.07) is there some evidence for a stable fee 
phase. These results are summarized in Figure 11. 

All micellar crystals obtained in this work are perfect 
single crystals displaying a high degree of order. Even 
in a periodic simulation box with four unit cells along a 
side, a single additional micelle can disrupt the result- 
ing lattice significantly. The number of micelles is con- 
trolled by changing the number of polymers in the box, 
as the system displays a remarkable stability in the av- 
erage aggregation number of the micelles that form. In 
the equilibrium lattice, micelles arc closely packed with 
a significant amount of overlap between the coronas of 
neighboring micelles. Polymers in the micelles are highly 
stretched across the liquid hydrophobic core through 
the solvent and qualitatively well described within the 
strongly stretched approximation [17, 30]. 



Stability of the bcc lattice near the disorder 
transition 

Our simulations results show a strong preference for 
bcc lattices near the disorder transition. A similar re- 
sult has been experimentally observed for PS-PI diblocks 
[10], where it has been attributed to the fact that near 
the disordered phase, micelle aggregation numbers are 
small. The phase diagram of f-star polymer systems 
shows that fee lattices are only stable for large number 
of arms / > 60, while bcc lattices are favored when the 
number of arms is small [36] . By considering polymeric 
micelles as f-star polymers, where the number of arms / is 
given by / ~ 2(N agg ) then bcc lattices are favored when 
the aggregation numbers are small (N agg ) < 30. This 
argument, however, hinges on two key assumptions: the 
dynamic nature of micelles does not play a significant role 
and that the hydrophilic blocks are sufficiently long. The 
first assumption is already somewhat problematic given 
the importance of polymer transfer found in this work. 
As for the second assumption, a criteria establishing how 
long hydrophilic blocks should be has been put forward 
in Ref. [8], where it is shown that if the size of the corona 
is L c and the core radius R c , bcc lattices are favored for 



L c /R c > 1.5. Our results for the A10B7A10 and AqBjA^ 
yield L c /R c ~ 0.9 and L c /R c ~ 0.7. It is therefore not 
possible to attribute the stability of the bcc lattices ob- 
served in our simulations as being a consequence of the 
small aggregation numbers (N agg ) < 30. 

It is tantalizing to interpret the stability of the bcc lat- 
tice in terms of the Alexander-McTague (AM) scenario 
[37], where it was argued that bcc should be generally 
expected to be the stable phase near a (weakly first or- 
der) disorder transition. Subsequent analysis however, 
showed that cubic lattices other than bcc cannot be ruled 
out near the disorder transition [38, 39]. In Ref. [38, 39] 
it is shown that the characteristic property of bcc lat- 
tices is that their free energy is closer to the disordered 
state, thus suggesting that bcc lattices follow an Ostwald 
step rule [40], where the solid phase that nucleates first 
is the one whose free energy is closest to the disordered 
(or fluid) state. In this case, the complete crystalliza- 
tion process would require an additional step, where af- 
ter the bcc crystallites are formed, they gradually evolve 
towards the stable thermodynamic phase. Certainly, it 
follows from our results that fee lattices are difficult to 
obtain by MD for the systems we simulate, but at least in 
one system (AqB-jAq at ksT /e — 1.07) the fee lattice has 
been obtained reproducibly, and did not proceed through 
an intermediate bcc step. In addition, for this very same 
system, closer to the disordered state (ksT/e = 1.1), no 
fee structure was found to be stable. Although not com- 
pletely conclusive, our results are more consistent with 
the bcc as being a stable thermodynamic phase near the 
disordered phase. 

There are serious limitations in identifying micelles as 
simple particles, because as the disordered phase is ap- 
proached micelle aggregation numbers decrease and poly- 
mers become essentially free. So the disordered phase is 
not a simple liquid as it is assumed by AM and all subse- 
quent work. Similar analysis in polymer melts [41], shows 
that in the vicinity of the spinodal, the only possible 
phase with cubic symmetry is bcc. Beyond the spinodal, 
other structures are favored [42]. Based on the previous 
discussion, we attribute the stability of the bcc phases 
observed in our simulations as a reflection of the admit- 
tedly non-rigorous statement that bcc phases are usually 
favored near the disordered transition. We defer to future 
work to establish this result within a rigorous framework, 
where all the nuances involved in micelle formation are 
properly taken into account. 

Implications for Pluronic systems 

The A10B7A10 and A20B14A20 polymers discussed in 
this paper provide coarse-grained descriptions of Pluronic 
F127. All simulations have been carried out in very good 
solvent conditions for the A beads and at total volume 
fractions of 15-20%. Experimental results in this region 
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of the phase diagram are surprisingly disparate. Simple 
cubic [31], bcc [32] and fee [33] have all been proposed 
as the structure in this region. Very recently, on the 
basis of new experimental data, the situation has been 
thoroughly reviewed by Li ct al. [34] (although for sig- 
nificantly higher temperatures), see also Ref. [35], but 
without clear conclusive results. Our theoretical analysis 
clearly favors the bcc lattice close to the disorder transi- 
tion. The comparison of our results with the experimen- 
tal F127 system is more accurate at low temperatures (at 
20-25 C), where the water can be considered as a good 
solvent for PEO, as discussed in Ref. [18]. 

Outlook 

We have shown that MD allows a detailed investiga- 
tion of both the dynamics as well as the thermodynamic 
equilibrium of micellar crystals. Many studies have been 
performed by modeling micelles as point particles, where 
the complex structure of the micelles is accommodated 
through refined two-body potentials, either derived an- 
alytically or empirically. While successful in many sit- 
uations, two-body potentials do not account for the dy- 
namic nature of micelles, which play a critical role in 
determining the phases of the system, particularly near 
the disorder transition. 

There are a few areas where further work needs to be 
performed. First, all simulations in this work have been 
performed in good solvent conditions. MD with implicit 
solvents of different quality are also of great interest, es- 
pecially to determine the phases of Pluronic systems over 
a wide range of temperatures [18]. Next, the largest 
micellar crystal formed in this work contains 2134 poly- 
mers (57,618 beads). We were unable to find any order in 
larger systems, even after running as many as 50 million 
time steps. It is possible that significantly longer sim- 
ulations may be required for larger systems. Also, the 
range of applicability of MD is restricted to a relatively 
narrow range near the disorder transition, as schemati- 
cally shown in Figure 11. Finally, the role of the solvent 
will need to be investigated in more detail, as it is found 
in the Rouse vs. Zimm dynamics for simple homopoly- 
mers [26] . Future studies will be necessary to completely 
clarify and expand all these issues. 

The relevance of our study goes beyond pure sys- 
tems. In Ref. [43] for example, Pluronic polymers have 
been used to template the growth of an inorganic phase 
of calcium phosphate, aimed at creating new polymer 
nanocomposites with lightweight /high strength proper- 
ties or that mimic the structure of real bone. An under- 
standing of the pure systems is clearly a prerequisite for 
accurate models of polymer nanocomposites. We hope 
to report more on this topic in the near future. 
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